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ABSTRACT 


his  paper  proposes  a rule  for  determining  when  to  start  collecting 
data  in  a queueing  simulation.  The  rule  is  designed  to  reduce  depen- 
dence between  the  empty  (queue)  and  idle  (servers)  initial  conditions  and 
the  collected  sample  record.  The  rule  is  an  outgrowth  of  earlier  work  by 
Fishman  and  Moore  (1978)  and  relies  on  a comparison  between  a priori 
information  on  the  activity  level  (traffic  intensity)  and  a corresponding 
sample  estimate  computed  during  the  course  of  simulation.  Experiments 
with  simulations  of  the  M/M/c  queue  with  c = 1,2,4  and  p = .7,. 8,. 9,. 95 
reveal  that  the  rule  reduces  and  in  most  cases  removes  the  dependence  on 
the  empty  and  idle  initial  conditions.  In  particular,  the  rule  begins 
data  collection  when  the  simulation  is  in  a congested  state  or  in  the 
steady  state.  The  rule  is  well  behaved  in  that  it  has  low  probabilities 
of  requiring  long  runs  before  data  collection  is  started.  Although  our 
data  suggests  an  association  between  the  rule's  performance  and  activity 
level,  the  performance  is  insensitive  to  variation  in  the  number  of  servers. 
Since  the  rule  is  based  upon  the  activity  level,  a parameter  that  frequently 
can  be  computed  from  the  input  parameters  of  the  simulation,  the  rule  is 

easily  generalized  to  a wider  class  of  queueing  simulations.-.  A subsequent 

\ 

study  (Adlakha  and  Fishman  1979)  demonstrates  the  appeal  of  the  starting 
rule  when  used  with  a proposed  stopping  rule  for  computing  interval  estimates 
of  parameters  of  interest. 


1.  Introduction 


■ 
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This  paper  proposes  a rule  for  starting  data  collection  in  a queue- 
ing simulation.  The  rule  is  an  outgrowth  of  work  by  Fishman  and  Moore 
(1978)  who  described  a framework  for  research  into  the  problem  of  reduc- 
ing the  effect  of  initial  conditions  on  sample  records  collected  for 
statistical  analysis  in  a discrete  event  simulation.  Initial  Conditions 
refer  to  the  states  of  critical  variables  at  the  beginning  of  a simula- 
tion run.  Because  of  the  dependence  among  phenomena  in  a simulation 
and  the  temporal  nature  of  much  of  this  dependence,  the  choice  of  ini- 
tial conditions  influences  the  observed  time  paths  of  these  phenomena. 

The  suggested  framework  in  the  Fishman  and  Moore  paper  encourages  the 
use  of  ancillary  information  such  as  the  theoretical  activity  level  in 
a queueing  simulation  to  provide  guidance  as  to  when  data  collection 
for  analysis  should  begin  in  a simulation.  In  particular,  if  a comparison 
between  the  theoretical  and  sample  activity  levels  in  a simulation  shows 
that  certain  conditions  are  met,  empirical  evidence  in  their  paper  shows 
that  the  queueing  simulation  is  in  a congested  state  for  an  important 
range  of  theoretical  activity  levels.  Therefore,  data  collection  initiated 
at  that  point  reflects  initial  conditions  for  a congested  system  rather 
than  for  an  undercongested  system,  as  would  occur  if  data  collection  starts 
at  the  beginning  of  a simulation  whose  initial  conditions  include  all 
idle  servers  and  empty  queues.  Since  most  simulators  are  inclined  to 
accept  upward  bias  due  to  initial  congestion  more  readily  than  downward 
bias  due  to  initial  undercongestion,  this  new-found  ability  to  induce 
congestion  is  a major  step  forward  in  alleviating  the  often  asserted 
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problem  of  initial  conditions:  When  does  a simulation  achieve  represen- 
tative state  behavior  of  the  system  being  studied? 

Ideally,  one  would  like  to  begin  data  collection  when  the  system 
is  in  the  steady  state.  However,  the  very  act  of  using  information  on 
the  sample  path  up  to  the  point  at  which  a decision  is  to  be  made  about 
starting  data  collection  from  that  point  onward  necessarily  creates  a 
conditionality  and,  therefore,  prevents  a determination  of  when  the 
steady  state  arises  in  a given  run.  Moreover,  from  a probabil ists 's 
point  of  view  one  regards  the  steady  state  as  a relative,  rather  than 
an  absolute,  concept  that  is  achieved  as  a limiting  operation.  In 
support  of  these  assertions,  the  empirical  evidence  in  the  Fishman  and 
Moore  paper  shows  that  the  steady  state  is  not  a point  of  attraction 
for  the  type  of  data-based  rules  they  employ.  However,  the  inducement 
of  congestion  is  possible. 

Fishman  and  Moore  presented  an  explicit  rule  for  achieving  this 

-f. 

congestion.  Although  the  rule  performed  well,  the  highly  variable  data 
collection  starting  times  to  which  it  led  allowed  for  the  possibility  of 
excessive  cost  before  achieving  congestion.  Their  report  concludes  with 
a recommendation  that  further  search  of  rules  to  induce  congestion  should 
consider  an  iterative  rule  which  they  describe.  This  paper  presents  the 
results  of  a considerably  more  elaborate  sampling  experiment  designed  to 
evaluate  this  rule,  hereafter  referred  to  as  Rule  1.  By  way  of  summary, 
the  results  to  be  presented  here  show  that  the  candidate  rule  can  induce 
congestion  with  considerably  less  risk  of  excessive  cost  for  a variety  of 
queueing  models  and  activity  levels. 

+The  rule  Is  called  Rule  2 in  their  paper. 
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Section  2 formulates  the  problem  and  Section  3 describes  the 
experimental  design  and  measures  of  evaluation.  Section  4 presents  an 
analysis  of  the  empirical  results  and  recommends  specific  values  for  the 
parameters  of  the  rule.  Section  5 presents  an  analysis  of  the  distri- 
bution of  starting  times  and  shows  that  the  new  recommended  rule  leads 
to  considerably  less  skewness  than  the  rule  recommended  in  Fishman  and 
Moore  (1978).  The  implications  for  cost  are  discussed  in  detail.  Section 
6 draws  overall  conclusions  and  describes  a proposed  algorithm  for  imple- 
menting the  rule. 

Before  beginning  the  technical  discussion,  we  wish  to  stress  two 
points.  Firstly,  the  ancillary  information  used  for  decision  making  here 
is  the  activity  level  of  the  queueing  system  being  simulated.  We  contend 
that  this  quantity  can  frequently  be  computed  from  the  input  parameters  of 
the  simulation.  This  computation  is  also  possible  in  many,  though  not  all, 
large  scale  queueing  simulations.  Although  we  recognize  that  our  results 
demonstrate  the  value  of  the  proposed  starting  rule  with  relatively  simple 
systems,  we  feel  that  its  applicability  with  perhaps  some  modification  to 
larger  scale  systems  is  a possibility  that  a thoughtful  simulation  user 
would  want  to  consider  seriously. 

The  second  point  concerns  the  ultimate  purpose  of  starting  data 
collection  in  the  steady  state  or  in  a congested  state.  In  a second 
paper  (Adlakha  and  Fishman  1979),  we  develop  a stopping  rule  for  data  col- 
lection which  when  used  with  the  starting  rule  proposed  here  produces 
confidence  intervals  via  the  autoregressive  method  (Fishman  1971)  whose 
coverage  rate  is  considerably  higher  than  that  reported  in  the  original 
Fishman  paper.  The  relatively  comprehensive  analysis  in  the  second  Adlakha 
and  Fishman  paper  makes  clear  the  attraction  of  using  the  suggested  start- 
ing and  stopping  rules  together. 
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2.  Problem  Formulation 

Consider  a simulation  model  of  a queueing  system  with  c servers 
in  parallel,  independent  interarrival  times  with  mean  1/X  and  indepen- 
dent service  times  with  mean  1/w  . Let  denote  the  elapsed  time 
between  arrivals  of  jobs  1-1  and  i and  let  denote  the  service 
time  of  arrival  i . Assume  that  the  simulation  begins  with  the  arrival 
of  job  1 to  an  empty  queue  and  c idle  servers.  Let  denote  the 
system  time  of  completion  i where  system  time  denotes  waiting  time  plus 
service  time.  Assume  that  an  ultimate  objective  of  analysis  is  to 
infer  the  characteristics  of  the  system  time  stochastic  process  from 
a sample  record  of  system  times.  Also  assume  that  given  the  choice 
between  starting  data  collection  in  an  undercongested  or  a congested 
system,  one  prefers  the  congested  one. 


Rvleetiny  a Rule 


After  n completions  occur  during  a simulation  run  one  can  esti- 

-1  n i n 

mate  1/X  and  l/u>  by  n } T.  and  n 1 £ S.  respectively. 

i=l  1 1*1  1 

These  estimates  are  unbiased  and  independent  of  initial  conditions. 

Now  the  activity  level  or  traffic  intensity  for  a queueing  system  usually 
is  defined  as 

p * arrival  rate/no.  of  servers  * service  rate 
= X/cw  , 


for  which  one  estimate  is 


= V,  s'/c  r,-,  T‘ 
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Since  pn  is  usually  a biased  estimator  of  p an  alternative,  presum- 
ably more  desirable,  estimator  is 


p P 

« - n 

•‘«v  ' 

since  E(?5n)  is  in  principle  derivable  for  most  common  interarrival  and 
service  time  distributions.  For  example,  in  the  case  of  exponential 
interarrival  and  service  times  E(pn)  = pn/(n-l)  so  that 


Pn  = (n-l)pn/n  . 


Let  S|  denote  the  condition  |pn  - p|  s <5  for  0 < 6 and  , 
the  condition  p^  > p ^ . Then  Fishman  and  Moore  recommend  that  data 
collection  begin  with  system  time  T+l  where 

Rule  ?.  T = min(n:  S1  and  S,,  hold). 

Condition  is  an  accuracy  requirement  designed  to  guarantee  that 
the  sample  activity  level  pn  is  acceptably  close  to  the  theoretical 
p . Condition  S?  is  a directional  criterion  designed  to  Insure 
that  data  collection  can  begin  omy  when  the  congestion  level  is  increas 
ing.  Intuitively,  obtains  if  job  n has  a large  service  time,  has 
a small  interarrival  time  or  has  both. 

Using  this  rule  on  an  M/M/1  queueing  model  Fishman  and  Moore 


conclude: 


On  the  basis  of  the  accumulated  empirical  evidence  to 
date,  one  inclines  to  recommend  the  use  of  rule  2 with 
5 * .0025  . Although  we  do  not  quarrel  with  this  recommen- 
dation, this  advice  should  be  regarded  as  a temporary  measure 
on  at  least  three  grounds.  Firstly,  we  have  no  experience 
with  p > .9  . Secondly,  we  have  no  experience  with  multi- 
server  systems.  Thirdly,  the  sample  quantiles  of  starting 
time  for  rule  2 and  6 * .9  In  ...  are  cause  for  concern. 
...  although  90  percent  of  the  starting  times  are  less  than 
3099  , one  percent  exceeds  48334  . In  our  opinion  the  risk 
of  excessive  cost  is  far  too  great  to  regard  rule  2 as  an 
end  In  itself. 


As  a candidate  for  Improved  performance,  they  suggest  beginning 
data  collection  with  system  time  T*1  where 

* <r  * 

Rule  1 T = min(n:  S^,  S?  and  S-j  hold). 


V l£I.n  ' pl  s s 

S2:  pI,n  * °I,n-l 
* 

S^:  ml  * n-1 


m * integer  > 0 , 


I * l(n-l)/mj 


In  particular  for  exponential  interarrival  and  service  times  one  has 


l s. 

j * ml+l  3 
n 

c ‘ l T, 

j=ml+l  J 


n - ml  - 1 
n - ml 


(1) 


i 
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In  words  this  rule  requires  one  to  use  a sample  activity  level  based  on 


at  most  m past  completions.  The  quantity  I denotes  the  number  of 


times  one  needs  to  reset  Pj  ^ ; i.e.,  the  number  of  iterations  minus  one. 

A little  thought  shows  that 


Pr(l  - D * 0 - %,>' 


i = 0,1 ... . 


where  q^  is  the  probability  of  success  on  a given  iteration.  Then  I 


has  a geometric  distribution  with  mean  (1  - qm)/qm  and  variance 


(1  - %,)/%,  • Also,  the  mean  number  of  completions  E(T  ) required  to 


meet  Rule  1 satisfies 


,(1  - %,>'%,  <E<n s m/%, 


Now  a user  may  choose  m to  suit  one's  convenience.  However,  from 


the  viewpoint  of  optimality,  one  prefers  the  m that  minimizes  m/q^  . 
If  several  m's  lead  to  identical  minima,  then  one  prefers  the  largest 


among  them  since  this  minimizes  var(ml)  = (1  - q^fm/q^  • In  the  next 


section  we  examine  the  performance  of  Rule  1 for  selected  values  of 


5 , m and  p on  simulations  of  the  M/M/1,  M/M/2  and  M/M/4  queueing 


systems. 


3.  E 


al  Design  a 


To  study  Rule  1,  we  used  a simulation  of  the  M/M/c  queueing  system 


with  the  experimental  design: 


p = .7,  .8,  .9,  .95 


c = 1,  2,  4 


6 = .001,  .0001 


m = 1000,  5000  . 


iSlj-','*  ' .V- . . r. - 


■ 


To  reiterate,  p denotes  the  activity  level;  c,  the  number  of  servers; 
6,  the  tolerance  to  be  achieved  before  comnenclng  dria  collection;  and 
m,  the  point  at  which  an  Iteration  terminates. 

For  all  simulation  runs  the  arrival  rate  was  fixed  at  unity  and 
the  mean  service  time  was  taken  as  c*p  . For  each  of  the  4*3*2x2  = 48 
experiments,  1000  independent  replications  were  performed.  Indepen- 
dence among  replications  was  achieved  by  using  the  last  seeds  of  random 
number  streams  for  a replication  as  the  initial  seeds  for  the  correspond 
ing  random  number  streams  in  the  next  replication.  On  each  replication 
of  each  design  point,  the  recorded  data  were: 

T = startinq  time  (number  of  observations  to  satisfy 
the  rule) 

I ■ number  of  iterations  minus  one 

*Tfl  system  t,me  (<lueue’n9  time  + service  time)  of  the 
(T+l)st  customer 

p,  j = activity  level  when  the  rule  is  satisfied  (computed 
using  (1)). 

Method  of  Eikiluation 

Let  X denote  a system  time  drawn  from  the  steady-state  distribu- 
tion of  system  time  for  the  M/M/c  queueing  system.  As  discussed  in  the 
Fishman  and  Moore  paper,  comes  from  a congested  system  if 

stochastically  dominates  X . One  says  that  the  random  variable  W 

with  cumulative  distribution  function  (c.d.f.)  Fu  stochastically  domi- 

w 

nates  the  random  variable  V with  c.d.f.  Fy  if  F^(u)  - Fy](u)  is 
non-negative  and  not  identically  zero  on  the  open  interval  0 < u < 1 
where  F^  and  F^  denote  the  right  continuous  inverses  of  F^  and 
Fy  respectively. 
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The  test  procedures  to  be  used  relate  to  the  hypotheses: 


Hq.-  Xyfl  comes  from  the  steady-state  distribution  of  system 
time  X . 

Hj : Xy+1  stochastically  dominates  X . 


Let  Tj  denote  the  starting  time  on  replication 
the  system  time  of  completion  T.+l  on  replication  j 

J 

replications  with  sample  data  X, 




dent  and  identically  distributed  when  using  Rule  1 . 
has  the  c.d.f. 


j a"d  Xj.Tj+l 

sJ 

. Consider  J 
which  are  indepen- 
Under  Hn:  X.  T 

j+1 


F ( x ) * 1 - e'(u,_A)x 


for  M/M/1 


F ( x ) = 1 ♦ 


(1  -P) 

TT+pITZp-IT 


-0)X 


2pc 

TTvpJTTFTj 


-(2u)-X)x 


for  M/M/ 2 


and 

F(x)  " 1 + ^ 4p^-P33)  e"U)X 


y --(4u)-x)x 
4p  - 3 6 


for  M/M/4 


where 


y = 


32pH 

3+9p+12p2+8p3 


Table  1 presents  the  corresponding  means,  variances  and  coeffi- 
cients of  variation.  Moreover, 


vj  * r'(xj.yi> 


has  the  uniform  distribution 


G(y)  s y 


0 s y < 1 


81.0000  1 361.0000 
32.4809  | 362.6139 
371.1839 

1.0000 

.977 


The  empirical  c.d.f.  of  the  Y.  s is 

J 


5 J '(0,y]<V 


0 < y < 1 


where 


if  0 < Y.  < y 

J 


otherwise. 


To  test  Hq  against  all  alternative  hypotheses,  one  examines  the  devia 


tions 


A(y)  = G,(y)  - y 


0 s y s 1 


1 1 


or  functions  ot  these  deviations,  lor  the  Kolmogorov -Smirnov  (KS) 
goodness-of-f it  test  one  uses  the  statistic 

0 » sup|A(y) | 

y 

for  which  critical  values  appear  in  Miller  (1953)  and  Owen  (196?). 

A second  test  procedure  uses  the  chi -squared  test  statistic 
? K 

X*"  * JK  }'  [A(1/K)  - A(  1 / K - 1/K)] 

1-1 

which  for  Integer  K and  large  J has  the  chi-squared  distribution 
with  K - 1 degrees  of  freedom  under  HQ  . A third  test  procedure 
uses  the  Anderson-Oarl inq  (AD)  test  statistic 

W?  ■ J / l[A(y)]2/y(l-y)}  dy 

0 

which  is  particularly  sensitive  to  departures  of  Gj(y)  from  G(y) 
in  the  tails  of  the  distribution.  Critical  values  of  the  test  statistic 
appear  in  Lewis  (1961). 

To  test  Hj  against  HQ  one  can  use  the  KS  test  statistic 

D“  « -inf  A(y) 

y 

Dwass  (1950)  describes  an  additional  helpful  measure  of  discrimination. 
The  statistic 

u ■ J 't— .0)(A<y))  ** 

gives  the  proportion  of  Gj(y)  that  lies  below  G(y)  - y . Under 
U has  the  uniform  distribution  on  (0,1)  . If  Hj  is  true  one  expects 
U to  be  close  to  unity. 


1? 


4 . Evaluation 

2 2 

Table  2 shows  the  D , x and  W statistics  for  HQ  against 
all  possible  alternatives.  The  data  show  that  a large  number  of  statis- 
tics are  significant  at  the  five  percent  level.  Before  we  draw  any 
inferences,  it  is  appropriate  to  consider  the  issue  of  multiplicity.  By 
multiplicity  we  mean  that  when  one  obtains  independent  observations  of 
several  test  statistics  using,  for  example,  the  five  percent  signifi- 
cance level,  then  one  would  expect  that  five  percent  of  the  statistics 
may  be  significant  when  HQ  is  true.  However,  the  number  of  signifi- 
cant statistics  in  Table  2 is  much  larger  than  one  would  expect  from 
random  variation  under  Hg,  and  is  evidence  that  the  occurrence  of 
these  significant  statistics  cannot  be  explained  satisfactorily  by  multiplicity. 

The  data  in  the  table  provide  substantial  evidence  to  reject  HQ 

at  each  design  point  for  p ^ .9  , Notice  that  for  p * .95  all  values 
2 2 

of  D , x and  W exceed  the  corresponding  critical  values  at  each 
design  point  with  the  exception  of  design  points  with  parameters  6 = .0001 
and  m = 5000  . This  shows  that  HQ  is  rejected  for  p = .95  and  each 
c , except  in  the  case  of  6 = .0001  and  m = 5000  , for  which  the 
general  pattern  tends  to  favor  Hg  . We  discuss  this  point  shortly. 

Recall  that  our  primary  objective  is  not  merely  to  detect  a departure 
from  the  steady-state  distribution  (reject  Hg)  . In  particular,  we  want 
to  check  for  a specific  type  of  departure  from  the  steady  state,  the 
stochastic  dominance  of  over  X . For  this  purpose,  we  use  the  D" 

and  U statistics.  The  statistic  U gives  the  proportion  of  the  empiri- 
cal c.d.f.  under  the  theoretical  c.d.f.  (prevalence  of  stochastic  dominance; 
and  the  statistic  D*  gives  the  magnitude  of  maximum  deviation.  Table  3 
presents  the  D"  and  U statistics  for  Hg  versus  Hj  . The  main 
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observations  are: 

(1)  For  p * .7,  .8  both  the  D*  and  the  U statistics 
support  the  hypothesis  of  stochastic  dominance  of  XT+^ 
over  X at  all  design  points,  with  the  exception  of  the 

U statistic  at  one  design  point  (p  ■ .7,  c = 4,  6 = .0001 
and  m * 1000)  . Also  observe  that  the  negative  deviation 
O'  is  reduced  as  6 decreases  or  as  c Increases. 

(2)  For  p = .9  the  D‘  statistics  support  the  hypothesis  of 
stochastic  dominance  at  all  design  points,  but  the  U 
statistics  generally  show  this  behavior  only  with  m = 5000  . 

(3)  For  p = .95  both  the  D"  and  the  U statistics  fail  to 
support  the  hypothesis  of  stochastic  dominance  at  all  design 
points. 

Although  Rule  1 fails  to  achieve  stochastic  dominance  for  p = .95  , 
the  statistics  in  Table  2 generally  favor  Hq  for  p = .95  with 

6 = .0001  and  m = 5000  . Here  we  consider  the  possibility  that  X sto- 
chastically dominates  Xy+^  , implying  that  Xy+^  comes  from  an  undercon- 
gested system.  Formally,  we  consider  for  p = .95  the  alternative  hypothesis 

X stochastically  dominates  Xy+j  . 

To  test  Hq  versus  H£  for  p = .95  , we  use  the  KS  statistic 

D+  = sup  A(y)  . 

y 
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The  results  appear  In  Table  4 . The  data  support  for  all 

design  points  with  the  exception  of  design  points  (p  ■ .95,  c • 1, 

6 =■  .0001,  m * 5000)  and  (p  - .95,  c ■ 2,  6-  .0001,  m ■ 5000)  . 

At  these  two  design  points  the  steady  state  Is  supported.  The  selec- 
tion of  Hq  over  H.,  In  these  two  cases  Is  reassuring,  for  If  HQ 
Is  more  credible  than  Hj  , we  prefer  that  HQ  also  be  more  credible 
than  . Achieving  the  steady  state  Is  acceptable,  because  this 
assures  that  the  Initial  conditions  of  the  simulation  no  longer  play 
a role. 

The  empirical  evidence  suggests  that  the  parameters  6 * .0001 
and  m * 5000  give  the  most  satisfactory  results  over  the  entire  range 


1 


» 


j 


Table  4 


Testing  Hq  versus  H2  for  p»  .95 


D+( .0381 )a  | 

c 

6 

m 

1 

2 

4 

.001 

1000 

.1375* 

.1479* 

.1122* 

5000 

.0960* 

.1078* 

.0848* 

.0001 

1000 

.1257* 

.1062* 

.1164* 

5000 

.0261 

.0269 

.0423* 

^Critical  value  at  the  5*  level. 
♦Significant  at  the  5*  level. 
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of  activity  level  p considered  In  this  study  [.7,  .95]  . For 
p < .9  , Rule  1 with  these  parameter  values  appears  to  Induce  sto- 
chastic dominance.  Though  the  test  statistics  for  p * .95  fail  to 
support  the  hypothesis  of  stochastic  dominance,  the  steady  state  seems 
to  be  achieved  for  c = 1,2  . As  stated  previously,  this  is  acceptable. 
Therefore,  we  recommend  the  values  of  6 * .0001  and  m * 5000  for 
use  with  Rule  1 . 

Since  data  collection  starts  with  the  system  time  Xy+j  , it 
Is  of  Interest  to  study  the  bias  and  the  dispersion  of  Xy+1 
Tables  5 through  7 present  the  sample  mean  Xj-+^  , the  sample 
variance  a (Xy+j)  and  the  samP1e  coefficient  of  variation  v(Xy+1) 
computed  over  the  1000  replications  of  each  deslqn  point  for  the  M/M/1, 
M/M/2  and  M/M/4  queueing  simulations.  The  data  show  that  Xy+j  is 
greater  than  the  corresponding  theoretical  mean  in  experiments  where 
Rule  1 achieves  stochastic  dominance,  and  Is  generally  less  than  the 
theoretical  mean  in  other  cases.  To  test  the  statistical  significance 
of  the  difference  between  the  theoretical  and  sample  mean  for  Rule  1 
with  our  recommended  parameters,  6 * .0001  and  m * 5000  , we  consider 
the  hypothesis 

h3:  ^Tfl  has  the  steady-state  mean 
against  the  alternative  hypothesis 

H^:  XT+1  does  not  have  the  steady-state  mean. 

The  test  statistics  appear  in  Table  8 . The  data  show  that  is 


Q > 


Sample  Mean,  Variance  and  Coefficient  of  Variation  for 
System  Time  XT., 


Activity  Level  p 

Statistic  6 m 0.7  I 0.8  I 0l9  I 0^95 


XTH  .001  1000  3.416  5.002  9.052  13.412 

5000  3.629  5.394  9.788  14.625 

.0001  1000  2.905  5.022  8.839  14.558 

5000  3.249  5.415  10.404  18.625 

2(Xt+1)  .001  1000  8.022  17.523  58.365  124.404 

5000  9.762  23.786  86.014  187.868 

.0001  1000  6.008  18.446  50.249  160.000 

5000  8.007  22.655  90.903  277.835 

^(XT>1)  .001  1000  .829  .837  .844  .832 

5000  .861  .904  .948  .937 

.0001  1000  .844  .855  .802  .869 

5000  .871  .879  .916  .895 


Theoretical  Quantities 

u 2.333  4 9 19 

o2  5.44  16  81  361 

v 1111 


Q > 
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Table  6 

Sample  Mean,  Variance  and  Coefficient  of  Variation  for 
System  Time 

M/M/2 


Activity  Level  p 

Statistic 

6 

m 

0.8 

0.9 

0.95 

*T*1 

.001 

1000 

3.453 

5.393 

9.235 

13.403 

5000 

3.489 

5.873 

10.221 

14.816 

.0001 

1000 

3.178 

4.865 

8.683 

15.631 

5000 

3.421 

5.708 

10.452 

18.698 

°2(XTa-l ) 

.001 

1000 

8.651 

18.971 

52.532 

138.434 

5000 

8.236 

23.631 

85.478 

212.971 

1 

1000 

6.692 

15.314 

49.372 

165.322 

5000 

8.244 

26.280 

86.144 

275.240 

V(XT>1 ) 

.001 

1000 

.852 

.808 

.785 

.878 

5000 

.823 

.832 

.905 

.985 

.0001 

1000 

.814 

.804 

.809 

.823 

5000 

.839 

.398 

.888 

.887 

— 

Theoretical  Quantities 

P 

2.7451 

I 4.4444 

9.4737 

19.4872  j 

2 

0 

6.4278 

17.2247 

82.4809 

362.6139 

V 

.924 

.934 

.959 

.977 

f 

i 

( 1 

' 


* . 


< D 
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Table  7 

Sample  Mean,  Variance  and  Coefficient  of  Variation  for 
System  Tine  XTH 
M/M/4 


Statistic 


.001  1000 

! 5000 


(XT+1)|  .001 


v(XT+1)  1.001 


.0001 


Theoretical  Quantities 


| Activity  Level  p | 

OKS 

■a 

0.95 

4.486 

6.309 

10.407 

16.105 

4.575 

6.439 

12.244 

17.411 

3.826 

6.050 

10.122 

16.712 

4.498 

7.570 

11.773 

18.141 

12.842 

23.458 

58.120 

150.670 

13.153 

24.898 

94.514 

225.918 

9.356 

24.461 

56.175 

176.487 

13.759 

40.273 

95.986 

237.959 

.799 

.768 

.733 

.762 

.792 

.775 

.794 

.863 

.799 

.818 

.740 

.795 

.825 

.838 

.832 

.850 

3.8002 

5.5857 

10.6898 

20.7370 

11.5072 

23.6341 

90.3110 

371.1839  ! 

.893 

.870 

.889 

.929 

r " — ■ — — 
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Table  8 

Student  t Statistics*  for  XT+1 
& ■ .0001  and  m ■ 5000 


P 

c 

.7 

.8 

.9 

.95 

1 

10.237* 

9.401* 

4.657* 

-.711 

2 

7.444* 

7.797* 

3.333* 

-1.504 

4 

5.959* 

9.888* 

3.496* 

-5.322* 

*The  entries  are  (XT+1  - p)/l600  / o(XT+1),  the  critical 
value  at  the  5%  level  Is  1.96  . 

Significant  at  the  5X  level . 


significantly  biased  upward  for  p s .9  . This  is  consistent  with  the 

i 

stochastic  dominance  established  for  p s .9  . For  p * .95  and  c = 4 
a significantly  downward-biased  condition  for  *T+1  is  indicated. 

The  data  in  Tables  5 through  7 show  that  for  each  design  point, 
v(XT+i)  is  less  than  its  theoretical  value.  This  observation  is  of  impor- 
tance because  it  indicates  the  probability  of  getting  an  from  the 

tail  distribution  is  not  increased  over  the  corresponding  steady-state 
probability  and,  therefore,  the  upward  bias  in  *T+1  is  not  the  result 
of  including  a disproportionate  number  of  extreme  waiting  times.  Figure 
1 presents  the  sample  cumulative  distribution  function  of  XT+1  and 
the  c.d.f.  of  X for  the  case  of  p ■ .9,  c ■ 1,  6 * .0001  and 
m ■ 5000  . The  figure  explains  the  fact  that  Rule  1 concentrates  the 
probability  of  system  time  away  from  the  tails. 


■ ■ 


V 
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Figure  1 

c.d.f.'s  of  XTfl  and  X 
P * .9,  c • 1,  6 • .0001,  m ■ 5000 

A basic  assumption  inherent  in  the  development  of  Rule  1 is  that 
the  correlation  between  X.^  and  T becomes  small  as  6 decreases. 
This  was  first  observed  in  Fishman  and  Moore  (1978)  and  was  one  of  the 
motivating  factors  in  considering  Rule  1 . In  principle,  one  would  like 
*T+1  an<*  f k*  independent.  Then  there  would  be  no  need  for  con- 
cern that  a small  T on  a particular  run  would  produce  an  undercongested 
system  for  starting  data  collection.  Table  9 presents  the  sample 
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correlation  between  and  T for  each  design  point.  Notice  that 

for  given  m , there  is  a tendency  for  the  correlation  to  diminish  as 
6 decreases,  albeit  exceptions  exist. 

Under  the  null  hypothesis:  corr(Xy+^,  T)  = 0 , the  sample  corre- 
lation coefficient  asymptotically  has  the  normal  distribution  with  mean, 
zero  and  variance  1/J  for  J independent  replications.  Let  us  con- 
centrate on  the  recommended  design  parameters  6 * .0001  and  m = 5000  . 
Significance  occurs  at  the  five  percent  level  at  p = .9  for  M/M/1 
and  M/M/4  and  at  p = .95  for  M/M/1  and  M/M/2  . Although  room  for 
improvement  exists  for  creating  a rule  that  makes  X^+^  and  T indepen- 
dent, the  relatively  small  magnitudes  of  the  significant  correlations 
encourages  us  to  recommend  Rule  1 at  present  with  minimal  concern. 

In  a study  involving  four  factors  p,  c,  6 and  m , it  is  of 
interest  to  see  how  performance  is  affected  by  different  levels  of  the 
factors.  One  way  to  measure  performance  with  regard  to  stochastic  domi- 
nance is  to  observe  how  D'  and  U change  with  the  alternative  levels. 

An  analysis  of  variance  (ANOVA)  enabled  us  to  investigate  these  questions. 
To  bring  the  data  in  closer  conformity  with  the  assumptions  of  ANOVA, 

0"  and  U were  transformed  to  standardized  normal  variates  under  the 
assumption  that  HQ  was  true.  For  example,  under  Hq  <J>-1(U)  is  a normal 
deviate  where  <t>’^  is  the  inverse  function  of  the  normal  distribution. 

ANOVA1 s were  performed  for  the  transformed  D"  and  U statistics 
separately.  Because  of  space  considerations,  we  report  the  most  relevant 
results  for  fixed  6=  .0001  and  m = 5000  only.*  They  show  that  for  the 

^Complete  details  of  the  ANOVA  studies  appear  in  Adlakha  (1979). 
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Table  9 

Sample  Correlation  Coefficients3  between  X.^  and  T 


c 

6 

m 

Activity  Level  p 

.7 

.8 

.9 

.95 

1 

.001 

1000 

.084* 

.072* 

.076* 

.282* 

5000 

.059 

.031 

.113* 

.211* 

.0001 

1000 

-.049 

-.035 

.025 

.114* 

5000 

.043 

.038 

.108* 

.129* 

2 

.001 

1000 

.044 

.062* 

.092* 

.275* 

.028 

.078* 

.241* 

.0001 

1000 

.044 

.027 

. 1 

.088* 

5000 

-.004 

.073* 

■ ! . 

.120* 

4 

.001 

1000 

CVJ 

o 

1 

.073* 

.176* 

5000 

.013 

.026 

.205* 

1000 

.031 

.041 

CO 

o 

o 

1 

.166* 

in 

■i 

.010 

-.036 

.063* 

.053 

aThe  critical  values  with  1000  observations  at  the  5%  and 
1%  significance  levels  are  .062  and  .081,  respectively. 

♦Significant  at  the  5%  level. 


range  .7  < p <■  .95  performance  is  relatively  insensitive  to  c . This 
is  a gratifying  result,  for  it  suggests  a generality  in  the  performance 
of  Rule  1 with  regard  to  varying  the  number  of  servers  in  a model.  How- 
ever, the  analyses  indicate  an  association  between  performance  and  p 
that  reinforces  the  visual  observations  in  Table  3 . Although  more  work 
is  needed  to  remove  this  association,  we  continue  to  recomnend  Rule  1 

for  p < .95  because  of  its  demonstrated  success.  Naturally,  a rule 

that  performs  independently  of  p is  the  ultimate  goal. 

5 . Distribution  of  Starting  Time 

A desirable  characteristic  of  a starting  rule  is  that  it  not  require 
excessive  amounts  of  computer  time.  The  computer  time  required  by  a 
starting  rule  is  essentially  composed  of  a program  set-up  time  plus  a 
running  time  component  that  is  generally  proportional  to  the  number  of 
observations  generated.  For  all  practical  applications  the  set-up  time 
is  insensitive  to  changes  in  starting  rules  when  compared  with  the  run- 
ning time  component.  Therefore,  it  is  of  interest  to  study  the  starting 
time  T associated  with  our  starting  rule,  since  it  is  this  variable 

that  determines  the  cost  (in  computer  time)  of  the  rule. 

We  are  interested  in  studying  the  variation  in  the  mean  starting 
time  in  response  to  changes  in  the  parameters  p,  c,  6 and  m . Adlakha 
(1979)  contains  these  results  for  the  complete  experimental  design.  Here 
we  focus  on  the  distribution  of  T for  Rule  1 with  the  recommended  param- 
eters <5  = .0001  and  m = 5000  . 


We  first  discuss  quantiles.  The  100*p  percent  quantile  of  the 
distribution  of  j is  min[n:  pr(y  s n)  * p]  . Table  10  presents 

the  sample  quantiles,  which  appear  to  be  relatively  insensitive  to 
o and  c . Although  the  95  percent  quantiles,  which  tend  to  be 
about  10,000  , are  comparable  to  those  observed  by  Fishman  and 

Moore's  starting  Rule  2 , the  higher  quantiles  and  the  maximum  value 
have  decreased  drastically.  This  Indicates  that  our  recommended  rule 
has  substantially  reduced  the  skewness  in  the  distribution  of  T . 

The  sample  mean,  standard  deviation,  and  coefficient  of  variation 
of  T also  appear  in  Table  10  . The  data  show  that  the  mean  starting 
time  varies  between  2837  and  3387  . These  mean  values  occur  around 
the  65  percent  quantile  as  compared  to  the  90  percent  quantile  for 
Rule  2 in  Fishman  and  Moore  (1978)  . The  coefficient  of  variation  of 
T is  approximately  equal  to  one  in  each  case.  These  observations  again 
suggest  that  the  starting  time  distribution  is  not  as  skewed  as  the  dis- 
tribution obtained  in  the  earlier  work. 

To  see  the  influence  of  this  reduction  in  skewness  on  the  cost  of 
a simulation  run,  we  consider  the  case  of  p s .9  with  the  cost  function 

C(T)  3 Cq  + CjT  + C2^[T*,«°)^  ’ 

for  example,  and  let  T*  = 30,000  . A cost  function  of  this  type  arises 
when  one  runs  out  of  computer  time  or  allocated  space  during  a simulation  run 
and  has  to  run  an  experiment  over  again.  One  can  also  conceive  of  such  a cost 
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Table  10 

Sample  Quantiles  of  Starting  Time  for  Rule  3 
<S  * .0001  . m • 5000 
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c 

1 

2 

1 

Np 

100»|K^ 

.7 

.8 

.9 

.95 

.7 

.8 

.9 

.95 

.7 

.8 

.9 

.95 

1 

46 

50 

67 

68 

47 

53 

60 

69 

55 

54 

66 

35 

2 

82 

74 

101 

96 

92 

83 

104 

133 

79 

100 

97 

78 

5 

154 

149 

192 

198 

169 

169 

184 

235 

168 

167 

176 

169 

10 

252 

279 

283 

318 

282 

262 

303 

357 

284 

303 

287 

328 

15 

336 

363 

377 

452 

377 

398 

396 

483 

397 

436 

418 

423 

20 

415 

469 

491 

556 

478 

477 

520 

603 

511 

563 

529 

527 

25 

523 

574 

635 

680 

568 

585 

627 

746 

631 

710 

691 

659 

30 

649 

684 

763 

828 

679 

744 

761 

882 

778 

818 

833 

803 

35 

771 

800 

877 

1004 

823 

868 

921 

1044 

906 

993 

1002 

966 

40 

916 

958 

1062 

1169 

985 

1016 

1112 

1301 

1064 

1172 

1194 

1118 

45 

1105 

1196 

1271 

1360 

1186 

1186 

1286 

1516 

1253 

1412 

1398 

1370 

50 

1297 

1366 

1509 

1576 

1437 

1464 

1537 

1783 

1563 

1737 

1766 

1621 

55 

1577 

1679 

1897 

1926 

1677 

1809 

1832 

2166 

1902 

2067 

2125 

2006 

60 

1999 

2068 

2329 

2375 

1930 

2238 

2153 

2694 

2357 

2109 

2631 

2468 

65 

2579 

2664 

2850 

2906 

2366 

2835 

2700 

3489 

2828 

2696 

3352 

2955 

70 

3418 

3384 

3823 

3722 

3257 

3613 

3553 

4609 

3807 

3450 

4271 

4103 

75 

4895 

4628 

5097 

5146 

4464 

4887 

4798 

5489 

5135 

4662 

5263 

5184 

80 

5428 

5382 

5554 

5687 

5352 

5406 

5563 

5977 

5627 

5448 

5596 

5579 

85 

5966 

5806 

6006 

6393 

5663 

5822 

6188 

6640 

6113 

6013 

6082 

6215 

90 

6884 

6683 

7560 

7816 

6540 

6752 

7202 

7982 

7213 

7150 

7170 

7173 

95 

10498 

10224 

10490 

10609 

10531 

10137 

10490 

10877 

10207 

10606 

10398 

10542 

98 

13534 

12013 

15306 

12825 

12418 

11581 

15306 

13896 

13985 

12524 

12566 

14193 

99 

15919 

15255 

18870 

15916 

15751 

13199 

17748 

16762 

15587 

15698 

14058 

18091 

min 

9 

2 

11 

6 

18 

9 

33 

17 

12 

7 

2 

8 

max 

27005 

24715 

34605 

25886 

27391 

26409 

31787 

25516 

25768 

21444 

21950 

27454 

T 

2906 

2850 

3188 

3174 

2839 

2877 

3057 

3387 

3062 

3035 

3145 

3149 

3561 

3332 

4043 

3627 

3442 

3218 

3682 

3696 

3496 

3407 

3337 

3656 

V 

1.23 

1.17 

1.27 

1.12 

1.20 

1.14 

1.12 

1.14 

1.06 

1.16 

*The  quantities  oT  and  vT  denote  the  sample  standard  deviation  and  the  coefficient  of 
variation  respectively. 
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function  when  computer  reliability  Is  low.  For  Rule  2 In  the  Fishman 
and  Moore  paper,  one  has 

E[C(T)]  “ cQ  + c1 E( T)  ♦ c2*.02 

and  for  Rule  1 

E[C(T)]  * c0  + ClE(T)  + c2«.001  . 

Since  the  mean  starting  time  obtained  Is  approximately  the  same  with 
Rule  2 and  the  recommended  Rule  1 , a significant  reduction  in  the 
expected  cost  Is  achieved  when  c2  Is  much  greater  than  c^  . Avoiding 
long  tall  starting  time  distributions  Is  clearly  a desirable  objective. 

6.  Conclusions  and  Proposed  Algorithm 

The  empirical  evidence  of  this  study  strongly  Indicates  that  the 
use  of  Rule  1 results  in  the  starting  of  data  collection  when  a system  is 
congested  (for  p s .90)  or  Is  at  least  in  the  steady  state  (for  p = .95) 
Although  a firm  theoretical  basis  for  this  dilution  of  the  Influence  of 
empty  and  idle  initial  conditions  remains  to  be  developed,  we  believe 
that  the  use  of  Rule  1 Is  a reasonable  recomnendatlon  for  a wider  class 
of  queueing  simulations  other  than  those  tested.  The  supposition  here 
is  that  one  can  compute  the  theoretical  activity  level  exactly.  The 
particular  form  that  the  activity  level  assumes  is  of  course  a function 
of  the  system  being  simulated. 
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There  are  several  ways  by  which  a user  can  implement  Rule  1 in  a 

simulation  program.  Essential  steps  are  provided  in  algorithm  START  . 

Algorithm  START 

Let  A = current  sample  activity  level 

B = old  sample  activity  level 
I = number  of  iterations  minus  one 
6 = tolerance  level  (given) 
m = iteration  length  (given) 
n = number  of  completions  . 

1.  Start  the  simulation  in  the  empty  and  Idle  state, 

n «-  0 , I «-  0 , B«-0. 

2.  Simulate  until  next  completion  and  n «-  n+1  . 

3.  Compute  A based  on  last  n - ml  completions. 

4.  If  B * 0 , |A  - p|  < 6 and  A > B go  to  step  7 . 

5.  If  n * m(I+l)  , B ^ A , go  to  step  2 . 

6.  I «-  1+1  , B «-  0 , go  to  step  2 . 

7.  Begin  data  collection  at  next  completion. 
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This  paper  proposes  a rule  for  determining  when  to  start  collecting 
data  in  a queueing  simulation.  The  rule  is  designed  to  reduce  dependence 
between  the  empty  (queue)  and  Idle  (servers)  initial  conditions  and  the 
collected  sample  record.  The  rule  is  an  outgrowth  of  earlier  work  by  Fishman 
and  Moore  (1978)  and  relies  on  a comparison  between  a priori  Information  on 
the  activity  level  (traffic  Intensity)  and  a corresponding  sample  estimate 
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computed  during  the  course  of  simulation.  Experiments  with  simulations  of 
the  M/M/c  queue  with  c ■ 1,2,4  and  p ■ .7,  .8,  .9,  .95  reveal  that  the 
rule  reduces  and  In  most  cases  removes  the  dependence  on  the  empty  and  Idle 
Initial  conditions.  In  particular,  the  rule  begins  data  collection  when  the 
simulation  Is  In  a congested  state  or  In  the  steady  state.  The  rule  Is  well 
behaved  In  that  It  has  low  probabilities  of  requiring  long  runs  before  data 
collection  Is  started.  Although  our  data  suggest  an  association  between 
the  rule's  performance  and  activity  level,  the  performance  is  Insensitive  to 
variation  In  the  number  of  servers.  Since  the  rule  Is  based  upon  the  activity 
level,  a parameter  that  frequently  can  be  computed  from  the  input  parameters 
of  the  simulation,  the  rule  Is  easily  generalized  to  a wider  class  of  queueing 
simulations.  A subsequent  study  (Adlakha  and  Fishman  1979)  demonstrates  the 
appeal  of  the  starting  rule  when  used  with  a proposed  stopping  rule  for  com- 
puting interval  estimates  of  parameters  of  Interest. 


UNCLASSIFIED 

MCUNlTT  CLASSIFICATION  OF  THIS  FAOSfOAa*  Oat*  Imm« 


